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ABSTRACT 


In this thesis several incompressible oscillatory flow and flutter problems 
were investigated. A previously developed unsteady panel code for single airfoil 
bending torsion flutter analysis was compared to Theodorsen's classical theory. 
The panel code agrees with Theodorsen's bending-torsion flutter analysis for 
natural frequency ratios (o) h /co a ) less than 1.2. In addition, a two airfoil unsteady 

panel code was modified for one degree of freedom flutter analysis. Code 
verification was accomplished by first comparing flat plate theory to the unsteady 
aerodynamic force and moment coefficients and then using the equation of 
motion to determine regions of instability. The possibility of active flutter control 
was investigated by positioning a small control airfoil in front of a neutrally stable 
reference airfoil. Results show that the flutter boundary may be changed through 
the scaling, placement, or oscillation of a second airfoil upstream. A comparison 
with pitch damping curves published by Loewy confirms that the code is capable 
of predicting wake-induced airfoil flutter. 
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I. INTRODUCTION 


A. GENERAL 

In this thesis, two unsteady panel codes UPOTFLUT described by Riester 
[Ref. 1] and USPOTF2C developed by Pang [Ref. 2] are modified and verified. 

The UPOTFLUT (Unsteady Potential Flow and Flutter) was developed by 
Teng [Ref. 3] for unsteady inviscid and incompressible flow over a single airfoil. 
The code is based upon the panel method by Hess and Smith [Ref. 4] that 
includes simple harmonic motion of the airfoil which is continuously shedding 
vortices into the wake. Riester [Ref. 1] extended the single airfoil code for 
bending-torsion flutter analysis. 

The USPOTF2F (Unsteady Potential Flow) code [Ref. 2] is a two airfoil code 
that permits steady and unsteady analysis of airfoil-wake interaction between 
two airfoils. A phase subroutine is added to convert the time dependent lift and 
moment histories into the frequency domain. 

Building upon the work of Riester [Ref. 1] which focused mainly on 
aerodynamic verification of UPOTFLUT code, the two dimensional bending- 
torsion flutter problem is examined. Theodorsen [Ref. 5, page 11] presents flat 
plate flutter speeds for comparison to UPOTFLUT computations. Specific cases 
are selected to compare the aerodynamic forces and solutions of the flutter 
determinant. 

The two airfoil USPOTF2 code’s unsteady aerodynamics are verified for pitch 
and plunge. Then the instability of one degree of freedom pitching oscillations is 
chosen to explore active flutter control. Pitch damping dependence upon 
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upstream airfoil size and position is fully explored. Emphasis is placed upon the 
wake effect and the code's ability to predict wake-induced flutter. 

Whenever possible, the theoretical work of Theodorsen [Ref. 5], Smilg [Ref. 
6], Loewy [Ref. 7] and other numerical examples are used for comparison to 
verify accuracy and trends. However, for many of the computations performed 
no other analytical theories and applicable experimental data exist. 

B. SCOPE 

Chapter II contains bending-torsion flutter verification using the single airfoil 
UPOTFLUT code and Theodorsen theory. Aerodynamic verification of force and 
moment coefficients using the two airfoil USPOTF2F code is accomplished in 
Chapter III. Chapter IV develops the equations of motion for steady-state 
oscillations in one degree of freedom. An in-ground effect numerical simulation is 
completed by oscillating airfoils out of phase. The positioning and size of the 
control airfoil for active flutter control are investigated in Chapter V. Application 
to rotary wings along with a comparison to wake-induced flutter theory is also 
discussed. 
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II. SINGLE AIRFOIL ANALYSIS 


The results discussed and presented in this chapter were produced from the 
UPOTFLUT code. This single airfoil analysis builds on work by Riester [Ref. 1]. 

A. INTRODUCTION TO UNSTEADY PANEL CODE THEORY 

Given a two-dimensional airfoil, the governing equation to solve the inviscid, 
incompressible flow over the airfoil is Laplace’s equation. Since the two- 
dimensional Laplace equation is linear, the principle of superposition may be 
applied. The four elementary flows: uniform flow, source flow, doublet flow and 
vortex flow are used to determine the flow around arbitrary bodies. For non¬ 
lifting bodies a uniform flow and source flow are required. For lifting bodies the 
vortex flow, must also be included to provide circulation. The boundary 
conditions which must be satisfied are flow tangency at the surface and the Kutta 
condition at the trailing edge. 

Steady-state panel methods as described by Hess and Smith [Ref. 4] may be 
extended to airfoils in simple harmonic motion. The continuous change in the lift 
implies a continually changing circulation about the airfoil. The Helmholtz vortex 
theorem requires that any change in the circulation around an airfoil must be 
matched by the appearance of an equal counter-vortex or starting vortex to 
achieve constant total circulation in the flow field. The starting vortices are 
assumed to be shed from the trailing edge after each time step and move with the 
local flow velocity. 

Riester [Ref. 1], page 23, showed that for single airfoil analysis 3 cycles is a 
sufficient wake length for amplitude and phase analysis when Kp = 1.0. The 
benchmark chosen for computations is 3 cycles, 65 time steps per cycle. A 
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NACA0007 airfoil geometry was chosen for comparison to flat plate theory. A 
thinner airfoil does not provide accurate results due to insufficient number of 
panels to adequately resolve the leading edge and vortex interaction between 
top and bottom surfaces. 

B. TWO DEGREE OF FREEDOM EQUATIONS OF MOTION AND 

FLUTTER DETERMINANT 

If a positive concentrated load, on a cantilever wing of uniform cross section 
is applied at the leading edge of the wing tip, the tip cross section will be 
displaced upwards and twist, thus increasing the angle of attack. If a similar 
concentrated load is applied at the trailing edge, the tip cross section will flex 
upwards and the angle of attack will decrease. A point exists between the 
leading edge and trailing edge where a concentrated load may be applied without 
twisting or causing the airfoil to rotate. This point is called the shear center or the 
flexural center of the airfoil. The elastic axis is defined as the locus of shear 
centers of the cross section, thus it is a natural reference line. 

Forces acting upon the airfoil are the elastic restoring forces, inertial forces 
and aerodynamic forces. Displacement variables for two degree of freedom, 
simple harmonic motion are pitch (a) and plunge (h): 

C a s torsional stiffne c . ! of wing 
C h s stiffness of wing in flexure 
M = mass of wing per unit span 
S 0 s static moment referred to elastic axis 
I a = moment of inertia about elastic axis 
The elastic restoring forces are hCj, and aC a . The inertial forces respectively in 
pitch and plunge are: 

a I a + h S a 
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h M + a S a 

Applying equilibrium by equating the inertial forces to the external aerodynamic 
and structural elastic forces and assuming zero structural damping yields: 

a I a + h S a = Moment - aC a 
h M + a S a = Lift - hC h 
The natural frequencies in pitch and plunge are: 

co h = 

(0 o = 

Substituting into the equations of motion: 

hM + aS„ + hMtoJ = Lift 
al a + h S a + aI o 0 )£ = Moment 

For simple harmonic motion: 

h = h o e i(0t a = a 0 e iwt 
h = -<B 2 h o e ,t0t a = -<n 2 a 0 e ia)t 
e 10 * (-to 2 h 0 M - to 2 a 0 S a + h 0 McoJ) = Lift 
e ,fflt (-co 2 a 0 I a -co 2 h 0 S a + a 0 I a £o 2 ) = Moment 

The panel code nondimensionalizes lift and moment using the characteristic 
length of chord (2b) vice semichord. 

Lift = (C Lo +Cu,)ipU 2 2b 

Moment = (C Ma + C Mh ) i pU 2 4 b 2 

Introducing dimensionless flutter parameters: 

Wing density parameter u. = ■ y 

Ttpb 
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Center of gravity from elastic axis: x a = —~ 

Mb 

Radius of gyration: r 2 = 

Mb 


Flutter frequency ratio: X=j 
The equations of motion become: 

[(-to*h„M -co 2 a 0 S a + h 0 Mo> 2 ) = (C^ + Cu,)± pU 2 2b] 
^2 [(-<o 2 «„I„ - to 2 h 0 S a +a„I<X) = (C Ma + Cmj, )|pU 2 4b 2 j 

Jk_ ax+ M^Y=fC +C 

b oa blcoj 1x1 u ' tcMcd 2 4b 2 
~ a<,r “ ~~b* a +a “ r “(lf) = ( CM “ +CMh * ,tM(0 2 4b 2 

-a 0 x 0+ |-^) -lj-tC^+Cu)^ 
a * r ‘((^) _1 )"^ Xo) = (Ckte +CMh) ^ 5 c 2 

Rearrange the two equations into the form: 

a(|) + B(a) = 0 

d(|) + E(a) = 0 


where 


l UJUJJ i) 


B = \ix a + 4C ^° 2 

jcaKp 
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D = (ix a + 


4C 


Mh 



E = \ir a 1- -*■ 


TtaKp 



A B 
|D E 



+ A 


Im aginary 


To find the critical flutter speed and frequency the characteristic equation is 
solved for U and to. Since A is complex both the real and imaginary parts must 
equal zero yielding two real quadratic equations for two unknowns. The real 
equation typically assumes a parabolic shape and the imaginary equation 



7 























C. MODIFICATION OF UPOTFLUT 

Significant modifications to the UPOTFLUT code involved 
nondimensionalizing flutter determinant parameters. 

The phase subroutine converts the time history of lift and moment into the 
form: 

F(t)=Amp* cos(a>t+<(>) 

Where Amp is the amplitude, (0 is frequency of oscillation and <j> is igle of 
phase lag of the aerodynamic forces to the motion of the airfoil. An iterative 
numerical scheme to find the phase shift involves summing the differences 
between actual time history and the fitted curve over a half cycle, then shifting <)> 

to minimize error. When iterations start with a phase error of n the total error is a 
maximum, hence the error curve’s derivative equals zero. If initial phase shift step 
size is too small, the numerical scheme will remain at a maximum instead of 
converging to minimize error. Initial step size of four degrees is sufficient to 
ensure convergence. 

D. VERIFICATION OF UPOTFLUT 

NACA TR-685 page 11 [Ref. 5] provides the closed form solution of a flat 
plate for comparison with the UPOTFLUT code. The results of Case (a) are 
identical to code validation performed by Riester [Ref. 1] page 101. Case (a) and 
Case (b) results are tabulated in Table 2.1. In general, the panel code results agree 
with Theodorsen’s theory. Next, the xalpha = 0.2 Case (h) was calculated using 
UPOTFLUT panel code. Case (h) physically represents a heavy wing with the 
center of mass aft of the elastic axis. The pivot axis is located at 35 percent of 
chord. These parameters make the cross section inherently susceptible to flutter. 
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As illustrated in Figure 2.3, the panel code results and Theodorsen theory diverge 
at natural frequency ratios greater than 0.8. 

The Theodorsen theoretical curve was verified using two dimensional 
flutter analysis described by Fung, page 235 [Ref. 8]. The analysis was 
completed by programming the fundamentals of flutter theory in MATLAB 
[Appendix A] using unsteady aerodynamic coefficients found on page 412 of 
Scanlan and Rosenbaum [Ref. 9]. The non-linear complex quadratic equations are 
solved by using a bisection root solving numerical scheme. Calculated non- 
dimensional flutter speeds coincided with published coefficients for Q> h /o) a less 
than 1.2. For co h /<o a greater than 1.2, the critical reduced frequencies 
encountered are greater than 4.0 based upon semichord (K t =4.0). At high 

reduced frequencies, the flutter determinant becomes ill conditioned and very 
sensitive to small changes in any of the flutter variables. Theodorsen did not 
have the computational power of computers in 1940 when TR-685 was 
published. Therefore, his results were recomputed. Although Theodorsen’s 
curve for Case (h) in TR-685, page 11 was found to be inaccurate for frequency 
ratios greater than 1.25 (as shown in Figure 2.3), the new theoretical curve using 
numerical methods did not coiTect the large discrepancy between his theory and 
the panel code. The source for this discrepancy may come from the computation 
of the unsteady aerodynamics or from the flutter calculations. Figures 2.4 
through 2.11 and Table 2.3 indicate close agreement between the aerodynamic 
coefficients for lower reduced frequencies with errors increasing slightly at higher 
frequencies. The panel code unsteady aerodynamic coefficients are derived and 
discussed by Riester [Ref. 1] and reproduced here specifically for TR-685 page 11 
Case (h). 
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The agreement between Theodorsen’s and the panel code unsteady 

aerodynamics suggests that the discrepancy in flutter coefficient is mainly due to 

the solution of the flutter determinant. The flutter coefficient is defined as: 

U 1 1 

bo)„ ” VX K, 

Figure’s 2.12 through 2.14 and Table 2.4 compare the real and imaginary 
roots of the flutter determinant from theory and the panel code. The intersection 
point of the real and imaginary roots represents the flutter frequency. It is shown 
in Figure 2.13 that the critical reduced frequency may differ by up to a factor of 
three from Theodorsen’s result for high critical frequencies and high frequency 
ratios. 
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AIRFOIL TYPE ! NACA 0007 AIRFOIL 
NLOWER - 100 , NUPPER • 100 


IFtAO NLOWER HUPPER 
0 100 100 
AIRFOIL TYPE 
3 

IRAHP IOSCIL ALPI ALPMAX PIVOT 

0 1 - 1.0 1.0 0.0 

FREQ RFOSTP RFQFWL 
0.04 0.02 0.2 

IGUST UGUST VGUST 

0 0 . 0 . 

ITRANS DELHX OELHY DELI PHASE 
0 0.00 .01 -.01 0.00 
CYCLE NTCYCLE TOL 

3 65 0.001 

naot 6 naot X aoa values Multiplied by 10 (integer) 

2 05 10 20 25 39 50 

wh/walpha ratio Maaa ratio xalpha (ralpha)**2 

0.4 3 0.2 0.25 

Consents... 

IRAMP 0: n/a RFREQ is based on full chord 

1: Straight ranp 
2: Modified ranp 

IOSCIL 0: n/a RFREQ is based on full cbord 

1: Sinusoidal pitch, notion starts at nin Aoa 

pivot : Position of elastic axis if leading edge-0 and trailing edge-1 

ITRANS 0: n/a 

1: Translational harnonic oscillation 

CYCLE : / of cycles for oscillatory notions 

-In case of ranp, cycle-1.5 denotes airfoil is held 
at nax aoa for the duration of .5 cycle 
-For steady state solution set it to 0 

NTCYCLE: # of tine steps for each cycle 

CYCLE*NTCYCLE is United to 200 currently. 

NAOT: # of input aoa for cp output 

- angles should be in Increasing order, 

- for oscilatory notions angles should increase 
first, then decrease. Decreasing angles are for 
the return cycle.. 


Figure 2.2 Typical UPOTFLUT.IN input file 
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TABLE 2.1 TR-685 CASE(A) and CASE(B) COMPARISON 
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Flutter Coefficient Comparison 
(TR-685 pg. 11 Case h Xalpha=0.2) 



Rgure 2.3 Rutter Coefficient Comparison 
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Real (Cl 


Re(Cl) Due to Pitch Comparison 
(TR - 685, Case h, pg. 11; 200 

panels; 3 cycles, 65 time steps 
per cycle) 



Kp = 2bw/U 


Figure 2.4 Real(Cl) Due to Pitch 
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Imaginary (Cl) 


Imaginary Part of Cl Pitch 
Comparison (TR -685, Case h, pg. 
11; 200 panels; 3 cycles 65 steps 

per cycle) 



Kp = 2bw/U 


Figure 2.5 Imaginary(Cl) Due to Pitch 
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Figure 2.6 Real(Cm) Due to Pitch 













Im(Cm) Due to Pitch Comparison 
(TR -685, Case h, pg. 11; 200 
panels; 3 cycles, 65 time steps 
per cycle) 


■o.oi 
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Re(CI) Due to Plunge Comparison 
(TR - 685, Case h, pg. 11; 200 
panels; 3 cycles, 65 time steps 
per cycle) 



Figure 2.8 Real(Cl) Due to Plunge 
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Figure 2.9 Imaginary(Cm) Due to Plunge 
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Figure 2.11 Imaginary(Cm) Due to Plunge 












TABLE 2.3 SINGLE AIRFOIL AERODYNAMICS 



































































































































































SQRT(X) 


Real and Imaginary Roots of 
SQRT(X) for wh/walpha = 0.4 



Figure 2.12 Roots for Frequency Ratio of 0.4 





















Real and Imaginary Roots of 
SQRT(X) for wh/walpha = 1.0 



Figure 2.13 Roots for Frequency Ratio of 1.0 
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SQRT(X) 


Real and Imaginary Roots of 
SQRT(X) for wh/walpha = 1.6 



Figure 2.14 Roots for Frequency Ratio of 1.6 
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TABLE 2.4 SOLUTIONS OF FLUTTER DETERMINANT 


The roots SQRT(X) of the real and imaginary equations against 1/Kt. 


TR • 685 pg. 11 Case h; mass ratio » 4; a - -0.3; Xalpha - 0.2;NACA0007 


itch: 0.01 h/2b plunge; 200 panels; 3 cycles 65 time steps 


wratio 0.4 


anel 


Reduced Fregue 



1/Ktheodorsen Real Panel 


0.25 


0.50 


0.83 


1.25 


1.67 


2.00 


2.50 


3.33 


5.00 


0.9628 


0.9771 


1.0273 


1.1422 


1.3145 


1.4963 


1.9345 



Imagina 


1.3573 


1.3533 


1.3386 


1.3315 


1.3338 


1.3400 


1.3536 


1.3689 


1.4982 


Real Theo 


0.9783 


1.0047 


1.0675 


1.1899 


1.3621 


1.5403 


1.9409 



Imagina 


1.3793 


1.3723 


1.3747 


1.3820 


1.3925 


1.4025 


1.4187 


1.4465 


1.4974 


wratio 1.0 


Reduced Fregue 


1/Ktheodorsen Real Panel 



0.25 


0.50 


0.83 


1.25 


1.67 


0.8235 


0.8293 


0.8577 


0.9223 


1.0026 


1.0602 


1.1184 


1.1628 


1.2075 


Imagina 


0.8843 


0.8743 


0.8642 


0.8607 


0.8631 


.8672 


0.8752 


0.8820 


0.9556 


theor 


Real Theory Imagina 


0.8317 0.8679 


0.8489 0.8661 


0.8880 0.8702 


0.9561 0.8772 


1.0303 0.8852 


1.0806 0.8917 


1.1293 0.9013 


1.1624 0.9160 


1.1767 0.9394 




I 

wratii 

o 1.6 


Reduced Fr 


1/Ktheodorsen 
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III. AERODYNAMIC COEFFICIENT VERIFICATION OF 

TWO AIRFOIL CODE 

The two airfoil USPOTF2F code was modified for one degree of freedom 
flutter analysis. This chapter verifies the aerodynamic force and moment 
coefficients by comparison to flat plate theory. 

A. INTRODUCTION TO TWO AIRFOIL PANEL CODE 

The USPOTF2F unsteady panel code is essentially the single airfoil UPOT 
code extended to two airfoils. Airfoil number one is treated as the master airfoil 
and remains at the origin of a cartesian coordinate system. Airfoil number two 
may be arbitrarily placed in reference to the master airfoil, as long as the leading 
airfoil’s wake does not impinge upon the trailing airfoil or its wake. All 
aerodynamic and flutter calculations pertain to the master airfoil. 

The phase subroutine computes the complex unsteady aerodynamic 
coefficients (refer to page 8), then the subroutine writes the coefficients to 
pitch.in and plunge.in files. Pitch.in contains K p , RealfC^), Imaginary 

Real(C Ma ) and Imaginary(C Mo ) tabulated in columns from left to right respectively 
for single degree of freedom pitch. Plunge.in contains the same information 
except the coefficients are due to single degree of freedom plunge in simple 
harmonic motion. The phase subroutine searches the last cycle of the time history 

of lift and moment for a maximum and minimum value. Average lift and moment 
coefficients are calculated from (C Lnux -C Lmin )/2. Each average is then 

subtracted from the time history respectively before curve fitting to a sine 
function: 

F(t)=Amp * sin(cot+<j>) 
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A function of this form assumes a zero average, or equal maximum and minimum 
magnitudes. Because, the pitch and plunge values start from zero at time equal to 
zero, thus this initial condition requires a sine function. The phase portion uses 
the first it radians of the last cycle for curve fitting. Selecting the positive area of 
the sine curve for integration avoids errors near k/2 and 3 jc/ 2 by not including 
positive and negative values in the integration interval. 

The product of cycles and time steps may not exceed 200 due to array sizing. 
As a warning to the user, note that the period of the last cycle is based upon the 
reduced frequency of the master airfoil. If the airfoils are oscillating at different 
reduced frequencies, the force and moment data for airfoil two in pitch.in and 
plunge.in files are invalid. 

To plot the wake vortices position and airfoil geometry, use the original 
USPOTF2C code. Fort. 10 and fort. 11 contain wake vortices position for each 
rime step. Cutting data from the last rime step in fort. 10 and fort.ll, and then 
pasting the data into separate files will allow for plotting applications. 

B. THEORETICAL COEFFICIENT ORIGIN 

The theoretical unsteady aerodynamic coefficient values for a flat plate were 
calculated using the Tables of Subsonic Incompressible Aerodynamic Coefficients 
from Scanlan and Rosenbaum [Ref. 9] pg. 412. Riester [Ref. 1 pg. 18] derives the 
complex coefficient of lift and moment values. The constants and equations were 
programmed in MATLAB [Appendix A] for comparison to USPOTF2F. 

C. GEOMETRY FOR SINGLE AIRFOIL COMPARISON 

Airfoils are separated by 50 chord lengths with no shift in the horizontal 
direction. The interaction between the two airfoils is assumed to be negligible at 
this distance. Inputs of 200 panels for each airfoil, pitch and plunge amplitudes of 
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one degree and 0.01 chord are used. Riester [Ref. 1] showed that small 
amplitudes yield more accurate results, especially in plunge at high reduced 
frequencies. As with UPOTFLUT aerodynamic verification, 3 cycles, 65 time 
steps per cycle was chosen as a benchmark wake length. 

D. AERODYNAMIC VERIFICATION 

The comparison between theory and USPOTF2F code is presented in Figures 
3.3 through 3.18 and Tables 2.1 through 2.4. For clarity, the computed 
aerodynamic coefficients are presented in phasor form as a magnitude and phase 
angle vice complex form. Phasor form allows for easier identification of errors. 

It is seen that the panel code results are in agreement with theory. 
Magnitude errors are less than five percent and phase error typically ranged from 
one to three degrees. A portion of the error may be attributed to comparison of 
an airfoil of finite thickness to a theoretical flat plate. 
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NUMBER OF LINES FOR TITLE 
1 

TWO HACA0007 OCSILLATIHG AIRFOILS 

I FLAG NLOWER NUPPER 

0 73 75 

NAIRFO, XSHIFT . YSHIFT, SCALE 
2 0 -100 1.0 

HACA AIRFOIL TYPE, 

7 

7 

ALFl ALP2 DALF1 DALF2 TCON1 TCON2 

0.0 0.0 0.0 0.0 0.0 0.0 

FREQ1 FRBQ2 FIVOT1 PIVOT2 

3.0 3.0 0.0 0.0 

UGUST VGUST DEL1IX1 OEL1IX2 OELilYl DELHY2 PHASE1 PHASE2 
0.0 0.0 0.0 0.0 0.1 0.1 0.0 0.0 
TF DTS1 DTS2 TOL TADJ SCL S CM SCAN NGIES 
6 65 0.0 .001 0.0 0.0 0.0 0.0 1 

STEADY OUTPUT 

false false 

STEADY—TRUE IF ONLY STEADY SOLUTION. FALSE OTHERWISE. 

OUTPUT—TRUE IF YOU WANT COMPLETE OUTPUT TO SCREEN. 

IFLAG : O IF AIRFOIL IS NACA XXXX OR 230XX 
1 OTHERWISE. 

HLOWER : NO. OF PANELS USED ON BOTH AIRFOIL LOWER SURFACES. 

NUPPER : NO. OF PANELS USED ON BOTH AIRFOIL UPPER SURFACES. 

NAIRFO t NUMBER OF AIRFOILS. 

XSHIFT : RELATIVE X DIST. FROM 2 AIRFOIL PIVOT POSITION WRT 
GLOBAL COORDINATE SYSTEM. 

YSHIFT ; RELATIVE Y DIST. FROM 2 AIRFOIL PIVOT POSITION WRT 
GLOBAL COORDINATED SYSTEM. 

HACA AIRFOIL TYPE : ENTER NACA 4 OR 5 DIGIT CODE FOR AIRFOILS., 

IF NOT A NACA AIRFOIL, SUPPLY AIRFOIL 
X(I),Y(I) COORDS. FOR BOTH AIRFOILS IN 
FILE CODE 2. 

ALF1/2 : INITIAL ANGLE OF ATTACK FOR AIRFOILS IN DEGREES. 

DALP1/2 : CHANGE IN AOA IN DEGREE FOR NON OSCILL. MOTION. 

MAX AMPLITUDE OF AOA IN DEGREE FOR ROT. HARMONIC MOTION. 

TCON1/2 t NON-DIMENSIONAL RISE TIME (Vinf.t/C) OF AOA FOR 
MODIFIED RAMP CHANCE IN AOA. 

FREQ1/2 s NON DIMENSIONAL OSCILL. (WC/Vlnf.) FOR HARMONIC MOTIONS. 
riVOTl/28 LENGTH FROM LEADING EDGE TO PIVOT POINT FOR LOCAL SYSTEM. 

(THE GLOBAL SYSTEM'S ORGIN IS THE FIRST AIRFOILS PIVOT POSITION 
UGUST : HAG. OF NON-DIM. GUST VELOCITY ALONG GLOBAL X DIRECTION. 

VGUST ! MAG. OF NON-DIM. GUST VELOCITY ALONG GLOBAL Y DIRECTION. 

DEMIX1/2: NON-DIM. TRANSLATIONAL CHORDWISE AMPLITUDE. 

DELIIY1/2: NON-DIM. TRANSLATIONAL TRANSVERSE AMPLITUDE (plunging). 
NIASE1/2! PHASE ANGLE IN DEGREE BETWEEN CHORDWISE AND TRANSVERSE 

TRANSLATIONAL OSCILL. WITH THE LATTER REF. TO THAT AIRFOIL. 

TF : FINAL NON-DIM. TIME TO TERMINATE UNSTEADY FLOW SOLUTION. 

DTS : STARTING TIME STEP FOR NON-OCIILL MOTIONS(TADJ-0). 

NO. OF COMPUTAIONAL STEPS PER CYCLE FOR HARMONIC MOTION 
(FOR 2 FREQ OCILL. IT USES THE LARGEST FREQ) 

BASELINE TIME STEP FOR ALL MOTIONS(TADJ NOT -0) 

0TS2 : STARTING NON-DIM TIME FOR SECOND AIRFOIL MOTION TO BEGIN. 

(O .0 BEGIN MOTION AT THE SAME TIME). 

TOL : TOLERANCE CRIERION FOR CONVERGENCE FOR (Uw)k and (Vw)k. 

TADJ : FACTOR BY WHICH DTS WILL BE ADJUSTED. 

SCLA : STEADY LIFT COEFF. FOR THE SINGLE AIRFOIL AT THE SPEC. AOA. 

SCM : STEADY MOMENT COEFF. FOR THE SINGLE AIRFOIL. 

SGAH : STEADY VORTICITY STRENGTH FOR THE SINGLE AIRFOIL. 

HGIES : OPTION TO CHANGE THE UNSTEADY KUTTA CONDITION. 

0 EQUAL PRESSURE AT THE TRAILING EDGE PANELS. 


Figure 3.1 USPOTF2C Input Namelist for Wake Analysis (Original) 
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NUMBER OF LINES FOR TITLE 
1 

TWO NACA0007 OCSILLATING AIRFOILS 

IF!AG HLOWER NUPPER 
0 100 100 
NAIRFO, XSHIFT , YSIIIFT, SCALE 

2 -43.98 -25 1.0 

UACA AIRFOIL TYPE, 

7 

7 

FREOl FREQ2 DALP1 DALP2 

0.1 0.1 1.0 1.0 

FREQSTEP1 FREQSTEP2 PIVOT1 PIVOT2 

0.0 0.0 0.0 0.0 

FIHALFREQ1 FINALFREQ2 DELHY1 DELHY2 

0.05 0.05 0.01 0.01 

CYCLES DTS1 0TS2 TOL TADJ SCL SCM SGAM NGIES 

3 65 0.0 0.001 0.0 0.0 0.0 0.0 1 

STEADY OUTPUT 

false false 

STEADY—TRUE IF ONLY STEADY SOLUTION. FALSE OTHERWISE. 

OUTPUT—TRUE IF YOU WANT COMPLETE OUTPUT TO SCREEN. 

IFiAG O IF AIRFOIL IS NACA XXXX OR 230XX 
1 OTHERWISE. 

NLOWER : NO. OF FANEI.S USED ON BOTH AIRFOIL LOWER SURFACES. 

NUPPER : NO. OF PANELS USED ON BOTH AIRFOIL UPPER SURFACES. 

NAIRFO : NUMBER OF AIRFOILS. 

XSHIFT ! RELATIVE X DIST. FROM 2 AIRFOIL PIVOT POSITION WRT 
GLOBAL COORDINATE SYSTEM. 

YSIIIFT ! RELATIVE Y DIST. FROM 2 AIRFOIL PIVOT POSITION WRT 
GIOBAL COORDINATED SYSTEM. 

NACA AIRFOIL TYPE s ENTER NACA 4 OR 5 DIGIT CODE FOR AIRFOILS., 

IF NOT A NACA AIRFOIL, SUPPLY AIRFOIL 
X(l),Y(I) COORDS. FOR BOTH AIRFOILS IN 
FILE CODE 2. 

DALP1/2 : CHANGE IN AOA IN DEGREE FOR NON OSCILL. MOTION. 

MAX AMPLITUDE OF AOA IN DEGREE FOR ROT. HARMONIC MOTION. 

FREQ1/2 : NON DIMENSIONAL OSCILL. (WC/Vin£.) FOR HARMONIC MOTIONS. 
PIVOT1/2: LENGTH FROM LEADING EDGE TO PIVOT POINT FOR LOCAL SYSTEM. 

(THE GLOBAL SYSTEM'S ORGIN IS THE FIRST AIRFOILS PIVOT POSITIOF 
DEUIY1/2: NON-DIM. TRANSLATIONAL TRANSVERSE AMPLITUDE (plunging). 

DTS : STARTING TIME STEP FOR NON-OCIILL MOTIONS(TADJ-0). 

NO. OF COMPUTAIONAL STEPS PER CYCLE FOR HARMONIC MOTION 
(FOR 2 FREQ OCILL. IT USES THE LARGEST FREQ) 

BASELINE TIME STEP FOR ALL MOTIONS(TADJ NOT -0) 

DTS2 : STARTING NON-DIM TIME FOR SECOND AIRFOIL MOTION TO BEGIN. 

(O TO BEGIN MOTION AT THE SAME TIME). 

TOL : TOLERANCE CRIERION FOR CONVERGENCE FOR (Uw)k and (Vw)k. 

TADJ : FACTOR BY WHICH DTS WILL BE ADJUSTED. 

SCLA : STEADY LIFT COEFF. FOR THE SINGLE AIRFOIL AT THE SPEC. AOA. 

SCM : STEADY MOMENT COEFF. FOR THE SINGLE AIRFOIL. 

SGAM ! STEADY VORTICITY STRENGTH FOR THE SINGLE AIRFOIL. 

NGIES : OPTION TO CHANGE THE UNSTEADY KUTTA CONDITION. 

0 EQUAL PRESSURE AT THE TRAILING EDGE PANELS. 

1 EQUAL TANGENTAL VELOCITIES AT THE TRAILING EDGE PANELS. 


Figure 3.2 USPOTF2F Input Namelist for Flutter Analysis (Modified) 
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Cl Magnitude (1.0 deg pitch; pivot 
about LE; NACA0007; 200 panels; 

3 cycles, 65 time steps per cycle) 



Figure 3.3 Cl Magnitude for Pitch about Leading Edge 
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Cl Phase (1 deg pitch; pivot about 
LE; NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 



Figure 3.4 Cl Phase for Pitch about Leading Edge 


















Cm Magnitude 


Cm Magnitude (1.0 deg pitch; pivot 
about LE; NACA0007; 200 panels; 

3 cycles, 65 time steps per cycle) 



Figure 3.5 Cm Magnitude for Pitch about Leading Edge 
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Cm Phase (rad) 


Cm Phase (1 deg pitch; pivot about 
LE; NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 



Figure 3.6 Cm Phase for Pitch about Leading Edge 
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TABLE 3.1 PITCH OSCILLATION PIVOT L.E. 







































































































































































































Cl Magnitude 


Cl Magnitude (1 deg pitch; 0.25c 
pivot; NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 



Figure 3.7 G Magnitude for Pitch about 1/4 Chord 















Cl Phase (1 deg pitch; 0.25 pivot; 
NACA0007; 200 panels; 3 cycles 
65 time steps per cycle) 



Figure 3.8 Cl Phase for Pitch about 1/4 Chord 















Cm Magnitude 


Cm Magnitude (1 deg pitch; .25c 
pivot; NACA0007; 200 panels; 

3 cycles 65 time steps per cycle) 



Figure 3.9 Cm Magnitude for Pitch about 1/4 Chord 
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PITCH OSCILLATION PIVOT 0.25 (1/4 cho 


bs; 65 time steps/cycle; 100 panels top and bottom 




theor 
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CL Phase 
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CM Phase 


anel 


3.81% 

•3.0664 


1.23% 

0.0055 
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3.13% 
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1.34% 

0.0083 

0.0076 1 

8.93% 

•1.4588 

-1.5436 

2.70% 

2.42% 

-2.7289 

•2.7809 

1.85% 

0.0111 

0.0102 ! 

i 8.36% 

-1.4219 

-1.4928 

2.26% 

1.64% 

-2.5638 

-2.6236 

1.90% 

0.0139 

0 

.0128 

7.63% 

•1.3854 

-1.4498 

2.05% 

0.98% 

-2.4135 

-2.4791 

2.09% 

0.0169 

0 

.0156 

7.95% 

-1.3495 

-1.4109 

1.95% 

0.15% 

-2.1614 

-2.2345 

2.32% 

0.0229 

0.0212 

7.52% 

-1.2793 

-1.3394 

1.91% 

1.70% 

-1.8061 

-1.8853 

2.52% 

0.0381 

0.0334 
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CL Real 


CM Real 
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•0.0522 
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0.0004 


0.0009 


0.0016 


.0026 


.0302 


0.0408 I 0.0222 
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0.0002 


0.0008 


0.0015 


0.0037 0.0025 
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TABLE 3.2 PITCH OSCILLATION PIVOT 1/4 CHORD 



























































































































































































































Cl Magnitude (1 deg pitch; 0.37c 
pivot; NACA0007; 200 panels; 

3 cycles, 65 time steps per cycle) 
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Figure 3.11 Cl Magnitude for Pitch about 0.37c 





















Cl Phase (rad) 



Figure 3.12 Cl Phase for Pitch about 0.37c 




















Cm Magnitude 


Cm Magnitude (1 deg pitch; 0.37c 
pivot; NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 



Figure 3.13 Cm Magnitude for Pitch about 0.37c 













Cm Phase (rad) 



Cm Phase (1.0 deg pitch; 0.37c 
pivot; NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 

o 

- 0.1 
- 0.2 
-0.3 
-0.4 
-0.5 
- 0.6 
-0.7 
- 0.8 
-0.9 
- 1 

0 2 4 6 8 

Kp = 2bw/U 

Figure 3.14 Cm Phase for Pitch about 0.37c 
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TABLE 3.3 PITCH OSCILLATION PIVOT 0.37C 


















































































































































































































Cl Magnitude 


Cl Magnitude (0.01 h/2b plunge; 
moment about LE; NACA0007; 
200 panels;3 cycles, 65 time 
steps per cycle) 



Figure 3.15 Cl Magnitude for Plunge 




















Cl Phase (rad) 


Cl Phase (.01 h/2b plunge; moment 
about LE; NACA0007; 200 panels; 

3 cycles, 65 time steps per cycle) 



Figure 3.16 Cl Phase for Plunge 



















Cm Magnitude 


Cm Magnitude (0.01 h/2b plunge; 
moment about LE; NACA0007; 200 
panels;3 cycles, 65 time steps 
per cycle) 



0 2 4 6 8 


Kp = 2bw/U 


Figure 3.17 Cm Magnitude for Plunge 













Cm Phase(.01 h/2b plunge; moment 
about LE; NACA0007; 200 panels; 

3 cycles, 65 time steps per cycle) 



Figure 3.18 Cm Phase for Plunge 
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TABLE 3.4 PLUNGE TRANSLATION 




















































































































































































































IV. SINGLE DEGREE OF FREEDOM FLUTTER 

When discussing flutter of a typical airfoil section, many text books begin the 
analysis of the flutter problem as a two-dimensional flexure-torsion problem. 
Smilg [Ref. 6] presents a theoretical study for one degree of freedom instability of 
pitching oscillations of an airfoil in incompressible flow. Theoretical results show 
that the single degree of freedom instability may occur if the airfoil’s moment of 
inertia is sufficiently large and if the pitch axis is located forward of the quarter 
chord. The pitch axis should not be too far forward of the leading edge so that 
the wing becomes stable. It is important to distinguish between stall flutter due to 
flow separation and inviscid incompressible attached flow flutter as discussed 
here. 

A. EQUATION OF MOTION FOR STEADY-STATE OSCILLATIONS IN 
ONE DEGREE OF FREEDOM. 

The classical spring-mass damper equation of motion is: 

mx + c x + kx = F 

The analogous equation for the pitch degree of freedom of a rigid airfoil in a two 
dimensional flow is: 

Ia + ga+ ka = M a 

Assuming no structural damping (g=0), k = Ico 2 , a = -(o 2 a 0 e ,a> \ and 
a = a 0 e ia *: 

I(co 2 -a> 2 ) = -M a 

Defining the aerodynamic pitching moment as M a = C ma ^-pU 2 4b 2 and 
substituting into the above equation gives: 
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For undamped oscillations of the airfoil to exist, the non-dimensional inertia 
parameter I 0 / 7 tpb 4 must exceed approximately 550. As the elastic stiffness 

increases, the inertia parameter must also increase for instability to occur. For a 
full scale aircraft at sea level, the inertia parameter typically will not exceed 30 for 
control surfaces, wings or stabilizers [Ref. 6]. 

For the oscillations to damp out, Imaginary(C Mot ) must be negative and the 

net work done on the airfoil per cycle must be negative. Negative work implies 
that the airfoil is performing work on the air, dissipating energy from the wing. 
When the aerodynamic moment and pitch displacement are in the same direction 
throughout the entire cycle positive work is performed on the airfoil. This 
transfer of energy is the origin of negative damping which leads to instability of 
pitch oscillations. The neutral stability boundary corresponds to 180 degrees 
phase lag between airfoil motion and lift and moment time history. This phase lag 
is defined by the necessary condition: 

Imag(C Mo ) = 0 

The real part of the lift and moment coefficients represents the in-phase 
aerodynamic loading while the imaginary part of the coefficients is the out-of¬ 
phase portion. 
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B. SINGLE AIRFOIL FLUTTER VERIFICATION 

Pitch damping (Imag(C Mo )) is plotted for several amplitudes in Figure 4.2. 

The critical reduced frequency representing the flutter boundary is approximately 
0.115 compared to a theoretical value of 0.078. The accuracy of the critical 
reduced frequency increases as the pitch amplitude is decreased. 

C. AIRFOIL IN-GROUND EFFECT SIMULATION 

To perform a ground effect simulation, the additional boundary condition of 
flow tangency at ground level requires a second airfoil’s reflection about ground 
level. Distance from the ground is non-dimensionalized by chord length. First, 
the steady state effects were examined by placing the two airfoils at zero degrees 
angle of attack and varying the vertical distance between them. At less than 3 
chord lengths of separation, the venturi effect or accelerated flow between the 
airfoils causes a decrease in lift corresponding to a positive pitching moment 
(Figure 4.3]. 

The critical reduced frequency of a pitching airfoil in ground effect initially 
increases slightly then decreases rapidly as the airfoil separation is reduced 
[Figure 4.4]. For an airfoil with a one foot chord, a natural frequency of one cycle 
per second and inertia parameter i/rcpb 4 equal to 550, the panel code calculates 

the flutter boundary to occur at a reduced frequency of Kp equal to 0.115 and 
C Mo equals -0.026613 out-of-ground effect. The out-of-ground effect flutter 
speed is 80 ft/sec and flutter frequency co/co a equals 1.464. At 4.5 chord lengths 

separation or 2.25 chord lengths above ground level the critical reduced 
frequency Kp equals 0.101 and C Ma equals -0.028942, which corresponds to an 
in-ground effect flutter speed of 125 ft/sec and flutter frequency (0/co a of 2.011. 
Clearly, the flutter speed and flutter frequency increase in-ground effect. 
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The same type of simulation was conducted for t\ o airfoils pitching in-phase. 
The flutter speed decreased and the critical reduced frequency of the flutter 
boundary increased [Figure 4.5]. 

A physical understanding of these results requires visualization of the wake. 
The in ground effect simulation tends to cancel the wake’s effect because the 
reflected wake has the same strength and opposite direction [Figure 4.6], In- 
phase oscillations will enhance the wake’s effect upon the aerodynamics of the 
airfoil [Figure 4.7]. The two wakes maintain the same separation as the airfoils 
have and the vortices have the same direction. 
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Figure 4.1 Phase Diagram for Airfoil Stability 








Cm Imaginary 


Im(Cm) for Various Pitch 
Magnitudes (pivot about LE; 
NACA0007; 200 panels; 3 
cycles, 65 time steps per cycle) 



Figure 4.2 Single Airfoil Flutter Verification 
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Steady State Flow of NACA 0007 
Airfoil in Ground Effect; 0 deg 
AOA; Moment about LE; 200 panels 



Figure 4.3 Steady State Aerodynamics, Moment Referenced to Leading Edge 
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Figure 4.4 Flutter Boundary In Ground Effect 
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TABLE 4.1 REFERENCES FIGURE 4.2 


Cm Imaainary due to Ditch 


NACA0007: 200 panels; 3 cycles 65 steps per cycle 


0.02 

0.000537 

0.001857 

0.003895 

0.04 

0.000639 

0.002247 

0.005146 

0.06 

0.000570 

0.002043 

0.004750 

0.08 

0.000393 

0.001456 

0.003459 

0.10 

0.000126 

0.000594 

0.001610 

0.12 

-0.000217 

-0.000480 

-0.000638 

0.14 

-0.000618 

-0.001718 

-0.003193 

0.16 

-0.001065 

-0.003084 

-0.005978 

0.18 

-0.001545 

-0.004544 

-0.008884 

0.20 

•0.002060 

-0.006101 

-0.012034 


TABLE 4.2 REFERENCES FIGURE 4.3 
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TABLE 4.3 REFERENCES FIGURE 4.4 



TABLE 4.4 REFERENCES FIGURE 4.5 
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ase Oscillations 



Figure 4.6 Out-of-Phase Oscillations, Wake Cancellation 


64 


Dimensional Time 

















Figure 4.7 In-Phase Oscillations, Wake Enhancement 


Wake Trace for In-Phase Oscillations (Kp = 0.1) 



Non-Dimensional Time 



















V. ACTIVE FLUTTER CONTROL 


The traditionally accepted methods of moving an airfoil’s flutter boundary is 
shifting the center of gravity or stiffening the structure. These two methods 
affect the inertial and elastic forces. In the past, aerodynamic forces were 
accepted as an uncontrollable quantity dependent upon reduced frequency. 
Through active flutter control, these forces can be controlled. 

An investigation of how the one-dimensional flutter boundary may be shifted 
by positioning an oscillating control airfoil upstream of a neutrally stable 
reference airfoil is presented in this chapter. The aerodynamic forces may be 
modified to stabilize an airfoil. 

A. INVESTIGATION OF AIRFOIL POSITIONING AND SIZE FOR 

ACTIVE FLUTTER CONTROL. 

The USPOTF2F code can accommodate and scale a second control airfoil 
while the reference airfoil remains fixed at the origin of the global coordinate 
system. The variables manipulated are horizontal spacing, vertical spacing and 
scaling of the contro * foil. The physical pitching frequencies of the two airfoils 
are equal for all cast presented in this chapter. The pitch axis is located at the 
leading edge. 

Pitch damping results are dependent upon the wavelength of the wake, 
defined as: 

v - 2nc 
Kp X 

where K p is reduced frequency based on full chord, c is chord length and X is 
the wavelength of the wake. Pitch damping refers to the imaginary part of the 
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moment coefficient. Positive Imag(Cm) corresponds to negative pitch damping or 
instability. 

Figure 5.1 displays the pitch damping-frequency dependence for identically 
sized airfoils staggered by two chord lengths ahead and below the reference 
airfoil. Figure 5.1 is very similar to Figure 4.2 constructed for single airfoil flutter 
verification. The staggered airfoil curve [Fig. 5.1] shifts down and slightly to the 
left when compared to Figure 4.2. 

Figure 5.2 shows the airfoil interaction at positive and negative Xshift values. 
For positive Xshift, the pitch damping asymptotically approaches the single airfoil 
value of 0.000126 [Table 4.1]. Between Xshift values of -5 to +10 the airfoil 
influence supersedes wake effects. The periodicity of pitch damping-Xshift is the 
wake wavelength. The mean of the sine curves in Figure 5.3 through Figure 5.5 
follows the single airfoil pitch damping curve [Fig. 4.2]. As the reduced 
frequency increases, the amplitude of Imag(Cm) decreases and the airfoil becomes 
more stable. In summary, the frequency dependence of pitch damping is inversely 
proportional to Cm imaginary period, mean and amplitude. 

The effect of scaling or sizing the control airfoil is displayed in Figures 5.6 
and 5.7. The ability of the control airfoil to affect the flutter boundary decreases 
proportionally with control airfoil size. As the size of an airfoil decreases, the 
strength of the wake core vortex shed from the trailing edge also decreases due 
reduced circulation around the control airfoil. 

B. APPLICATION TO UNSTEADY AERODYNAMICS OF ROTARY 
WINGS. 

The USPOTF2F code is well suited for unsteady aerodynamic simulation of a 
hovering helicopter. A two dimensional approximation considers the influence of 
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shed vortices that have translated below the rotor disc and passed under the next 
blade. Wake-excited flutter has been documented in laboratory tests and during 
initial whirl tower testing of the AH-56 Cheyenne. Wake-excited flutter will only 
occur at low inflow rates when wake spacing is minimal. The low inflow 
configuration occurs during ground operations when blade pitch angles are small. 
In an actual helicopter rotor flow field, the shed vorticity sheets will superimpose 
upon each other becoming stronger with the passage of each successive blade. 
A drawback to the USPOTF2F code is that it only accounts for the wake of one 
previous blade. 

C. COMPARISON TO TWO DIMENSIONAL APPROXIMATION TO 

WAKE FLUTTER OF ROTARY WINGS BY ROBERT LOEWY. 

1. Parameters 

In Reference 7, page 90, Figure 15 Loewy presents a graph of pitch 
damping coefficient versus frequency ratio as a function of inflow parameter [Fig. 
5.8]. The inflow parameter is the vertical spacing of the wakes of previous blades 
corresponding to Yshift in USPOTF2F code. The frequency ratio m=to/Q is 
divided into two portions, the integer and non-integer part. 

The non-integer portion of m represents wake phasing. When m equals 
0.5, this corresponds to the wakes being out of phase as shown in Figure 4.6. 
The vortex from the control airfoil or preceding blade is directly below a core 
vortex from the reference airfoil in the opposite direction. The wakes are counter 
phase or 180 degrees out of phase. When m equals zero or one, the phase shift is 
0 and 360 degrees respectively [Fig. 4.7]. The vortex wake sheet has equal 
spacing vertically with equal vortex strength and direction for identical airfoils. 
Couch in Reference 10, page 39 treats m exclusively as a phasing parameter in his 
finite wake theory. 
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The integer portion of m represents the reduced frequency of the airfoil. 


(0 



Multiplying by the inverse of aspect ratio for a rectangular wing yields reduced 
frequency. 


0 ) c cue c c 



Reduced frequency is physically defined as the ratio of chord length to wake 
wavelength as shown above. 

2. Comparison to Loewy Results 

Figure 5.8 shows areas of negative pitch damping between m=0.5 to 1.0 
and m=1.5 to 2.0. The m=0.5 to 1.0 region corresponding to a phase shift of 180 
to 360 degrees is illustrated in Figure 5.3 when Xshift ranges from -31.440 -62.8, 
which is a region of instability. Figure 5.3 through Figure 5.5 confirms that as the 
reduced frequency increases (Kp=0.1 to 1.0), the amplitude of the pitch damping 
curve will decrease. The difference in shape of Loewy’s pitch damping curve 
[Fig. 5.8] and sinusoidal shape of Figure 5.3 is due to Loewy’s use of an infinite 
number of wakes in his classical closed form solution while USPOTF2F code only 
accounts for one wake. 

D. USPOTF2F CODE LIMITATIONS 

Theodorsen’s flutter theory for simple harmonic motion assumes sinusoidally 
varying lift and moment coefficients, which are represented mathematically by 
e ,a *. The primary limitation is that the unequal physical frequencies of control 
and reference airfoil yield higher harmonic time histories of sine and cosine. The 
phase subroutine attempts to curve fit a sine wave to a higher harmonic time 
history resulting in rapid inconsistent phase shifts. 
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Cm Imaginary 


Cm Imaginary Values for 
Staggered Airfoils 
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Figure 5.1 Frequency Dependence of Pitch Damping Xshift=-2, Yshift=-2 
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TABLE 5.1 DATA REFERENCES FIGURE 5.1 


Pitching Airfoils in Phase at Various Freauencies 


Pivot about LE; NACA0007; 200 panels; 3 cycles 65 ste 


Xshift -2.0; Yshift -2.0; Oaipl - 1.0; Dalp2 - 1.0; 


Airfoils 1 and 2 freauencv sweep from .025 to .250; 




Freal 


0.025 


0.038 


0.050 


0.063 


0.075 


0.100 


0.113 


0.125 


0.138 


0.150 


0.163 


0.175 


0.188 


0.200 


0.213 


0.225 


0.238 


0.250 


Frea2 


0.025 


0.038 


0.050 


0.063 


0.075 


0.088 


0.100 


0.113 


0.125 


0.138 


0.150 


0.163 


0.175 


0.188 


0.200 


0.213 


0.225 


0.238 


0.250 


m 


Cm Imaaina 


0.000959 


0.000972 


0.000822 


0.000545 


0.000222 




-0.000549 


-0.000982 


-0.001432 


-0.001887 


-0.002345 


•0.00282 


.003297 


-0.003774 


•0.004234 


-0.004698 


-0.005149 


•0.005603 


-0.00607 
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Imaginary Cm 


Cm Imaginary Plot for Constant Kp 

0.1 

0.002 
0.0015 
0.001 
0.0005 
0 

-0.0005 
- 0.001 
-0.0015 
- 0.002 
-0.0025 
-0.003 

- 1 0 0 1 0 20 
Xshift (chord lengths) 

Figure 5.2 Pitch Damping at Small Xshift Values Compared to Wake Wavelength 
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TABLE 5.2 DATA REFERENCES FIGURE 5.2 










































































Imaginary Cm 


Cm Imaginary Plot for Constant Kp = 0.1 



Figure 5.3 Pitch Damping-Xshift Dependence 
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Imaginary Cm 


Cm Imaginary for Constant Kp = 0.2 



Figure 5.4 Pitch Damping-Xshift Dependence 
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TABLE 5.3 DATA REFERENCES FIGURE 5.3 AND FIGURE 5.4 


Pitching Airfoils 


Pivot about LE; NACA0007: 200 panels; 3 cycles 65 steps per 


Yshift -2.0: Dalpl 1.0: Dalp2 1.0 


Xshift varied from 0.0 to -120 chord lengths 


- 0.1 


Xshift 


-6.28 


-12.57 


-18.85 


-25.13 


-31.42 


-50.26 


-56.55 


-62.83 


-69.11 


-75.40 


-81.68 


-87.96 


-94.25 


100.53 


106.81 


113.10 


119.38 


Cm Imagina 


0.000737 


-0.001992 


-0.002929 


•0.002595 


- 0.001021 


0.001010 


0.002836 


0.003599 


0.003089 


0.001560 


-0.000474 


-0.002236 


-0.003109 


-0.002814 


-0.001184 


C.000944 


0.002751 


0.003464 


0.002851 


0.001307 


ni 


K 


Xshift 


0.00 


-3.14 


-6.28 


-9.42 


-12.57 


-15.71 


-18.85 


-21.99 


-25.13 


-28.27 


-31.42 


-34.56 


-37.70 


-40.84 


-43.98 


-47.12 


-50.27 


-53.41 


-56.55 


-59.69 


Cm Imagina 


-0.001999 


-0.005072 


-0.006196 


-0.005637 


-0.003424 


-0.000595 




0.002647 


0.001769 


-0.000400 


-0.003206 


-0.005419 


-0.006422 


-0.005997 


-0.003724 


-0.000667 


0.001671 


0.002466 


0.001482 


-0.000611 
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TABLE 5.4 DATA REFERENCES FIGURE 5.5 


Pitchina Airfoils Ko - 1.0 


Pivot about LE; NACA0007; 200 panels; 


Yshift *2.0; Dalpl 1.0; Dalp2 1.0 


Xshift varied from 0.0 to -18.22 chord lengths 


Xshift 


0.00 


-0.63 


-1.26 


- 1.88 


-2.51 


-3.14 


-3.77 


-4.40 


-5.03 


-5.65 


-6.28 


-6.91 


-7.54 


-8.17 


•8.80 


-9.42 


-10.05 


- 10.68 


-11.31 


-11.94 


-12.57 


-13.19 


-13.82 


-14.45 


-15.08 


-15.71 


-16.34 


-16.96 


-17.59 


-18.22 


5 cycles 39 ste 


-0.027219 


-0.0280 60 


-0.029022 


-0.029452 


.029082 


-0.027958 


•0.026568 


-0.025448 


-0.025049 


-0.025478 


•0.026566 


-0.027912 


•0.029073 


-0.029587 


•0.029219 


umnm 


-0.026656 


-0.025488 


•0.025106 


-0.025538 


-0.026627 


-0.027934 


-0.029100 


-0.029640 


-0.029299 


-0.028121 


.026633 


-0.025452 


.025088 


-0.025537 


3 cycles 65 ste 


-0.027008 


-0.027836 


-0.028742 


•0.029159 


•0.028802 


-0.027730 


-0.026411 


-0.025333 


-0.024928 


-0.025330 


-0.026357 


-0.027685 


-0.028777 


-0.029245 


•0.028880 


-0.027792 


mmn 


-0.025244 


-0.024793 


-0.025154 


m 





















































































































TABLE 5.5 DATA REFERENCES FIGURE 5.6 


Pitching Airfoils Kp - 0.1 


Pivot about LE; NACA0007; 200 panels; 3 cycles 65 steps per cycle; 


Yshift -2.0; Daipl = 1.0: Da!o2 « 1.0: 


Xshift varied from 0 to 120 chord lengths for scales of 0.25. 0.50 and 1 


ml 


Xshift 


0.00 


•6.28 


-12.57 


-18.85 


-25.13 


-31.42 


7.70 


-43.98 


-50.26 


-56.55 


-62.83 


-69.11 


-75.40 


-81.68 


-87.96 


-94.25 


100.53 


106.81 


113.10 


119.38 


Scale 1.0 


0.000737 


•0.001992 


-0.002929 


-0.002595 


- 0.001021 


0.001010 


0.002836 




0.003089 


0.001560 


.000474 


-0.002236 


-0.003123 


-0.002814 


-0.001184 


0.000944 


0.002751 


0.003464 


0.002851 


0.001307 


Scale 0.5 


0.000429 


-0.000939 


-0.001471 


-0.00125 


-0.000415 


0.000659 




0.002005 


0.001742 


0.000901 


-0.000172 


-0.001098 


-0.001599 


-0.001355 


-0.000495 


0.000597 


0.001536 


0.001932 


0.001641 


0.000765 


Scale 0.25 


0.000317 


-0.000383 


-0.000638 


-0.000532 


-0.000114 


0.000437 


0.000903 


0.001116 


0.000988 


0.000557 


0.000013 


-0.000455 


-0.000697 


-0.000581 


-0.000154 


0.000402 


0.000869 


0.001076 


0.000929 


0.000477 
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Figure 5.7 Yshift Dependence 
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TABLE 5.6 DATA REFERENCES FIGURE 5.7 
































































Figure 5.8 Loewy Pitch Damping Curves [Ref. 7] 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. SINGLE AIRFOIL ANALYSIS 

Aerodynamic verification of UPOTFLUT code is consistent with results by 
Riester [Ref. 1]. Incompressible bending-torsion flutter computations compare 
favorably with Theodorsen theory when the natural frequency ratios are less than 
1.2. It was shown that the critical reduced frequency may differ by up to a factor 
of three from Theodorsen’s result for high critical frequencies and high natural 
frequency ratios. 

B. TWO AIRFOIL ANALYSIS 

The nonlinear theory for simple harmonic motion, magnitude and phase 
relationships that exist between airfoil motion and aerodynamic forces have been 
verified by comparison to the classical Theodorsen analysis. Agreement with 
theory consistent with previous panel codes was obtained. Aerodynamic 
coefficient magnitude errors were less than five percent and phase errors typically 
ranged from one to three degrees. 

A single degree of freedom flutter analysis of an airfoil in-ground effect 
simulation showed that the flutter speed increases as ground distance decreases 
due to wake cancellation. The one dimensional flutter boundary may be shifted 
by positioning and oscillating a control airfoil upstream of a neutrally stable 
reference airfoil. An investigation of positioning and size of a control airfoil for 
active flutter control showed the pitch damping was periodic with wake phasing. 
Pitch damping period, mean and amplitude are inversely proportional to reduced 
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frequency. A trend comparison with the wake flutter theory by Loewy illustrated 
consistent code results. 

Theodorsen flutter theory for simple harmonic motion assumed sinusoidally 
varying lift and moment coefficients, which were represented mathematically by 
e l<1 *. The physical frequency differences between the control and reference 
airfoil yielded higher harmonics in the time history, causing Theodorsen theory to 
be invalid. 

C. RECOMMENDATIONS 

A drawback to this code was the lengthy computer run times which were 
approximately 60 hours per figure on a 100 Mhz SGI Workstation. Run times 
could be reduced by running the code in parallel. Dividing the simulation into 
separate directories or cases and running up to ten workstations simultaneously 
make time requirements manageable. This requires the availability of batch 
queues and access to a network. Streamlining of the code by an experienced 
programmer is recommended. 

The frequency dependence of the control airfoil upon a reference airfoil was 
not resolved. The non-sinusoidal lift and moment time histo. . will require 
computations to be completed in the time domain. The concept of active flutter 
control is feasible and certainly worthy of further theoretical study. 
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APPENDIX 


This appendix contains theoretical Theodorsen flutter analysis as described 
by Fung [Ref. 8]. The MATLAB program was written for the Macintosh. Next, 
time histories of lift, moment and pitch angle for various cases are illustrated. 
Finally, three sample wake position plots show the core vortex position at the 
final time step. 
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% Hus program performs theoretical theockirseh analysis as 

% described hy Fi^. 

clear 

clg 

% 

% Aerocfynamic Coefficients 
K(l) = 100; 

Ui(l) = 1; 

Maipha(l) = 0.375; 

Laipta(l) = 0.5; 

K(2) = 4; 

Ui(2) = 0.9848-0.2519*1; 

Maipiia(2) = 0.375-0.250*i; 

Lalpha(2) = 0.42179-0.49423*1; 

K(3) = 2; 

Lh(3) = 0.9423-0.5129*1; 
tolpl0(3) = 0.375-0.5*1; 

Lal|ilia(3) = 0.18580-0.98405*1; 

K(4) = 1.2; 

Ui(4) = 0.8538-0.8833*i; 

Malpha(4) = 0.375-0.83333*1; 

Ialptia(4) = -0.38230-1.59487*1; 

K(5) = 0.8; 

Ui(5) = 0.7088-1.3853*1; 

MhlplH(5J a 0.375-1.25*1; 

LalptH(5) = -1.5228-2.27119*1; 

K(6) = 0.6; 

Ui(6) = 0.5407-1.9293*1; 
telpha(6) a 0.375-1.66667*1; 
talpta(6) a -3.17490-2.83045*1; 

K(7) = 0.5; 

Ui(7) = 0.3972-2.3916*1; 

*felpha(7) = 0.375-2.0*1; 

Lalpha(7) = -4.8860-3.1860*1; 

K(8) = 0.4; 

Ui(Q) = 0.1752-3.1250*1; 

M?lrJ«(8) = 0.375-2.5*1; 

Lalpha(8) = -8.1375-3.5625*1; 

K(9) = 0.34; 

Ui(9) = -0.022-3.8053*1; 
tolplia(9) = 0.375-2.94118*1; 

Lalpha(9) = -11.714-3.7396*1; 

K(10) = 0.3; 

th(10) = -0.1950-4.4333*1; 
telplBUO) = 0.375-3.333333*1; 
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LalttaUO) = -15.473-3.7822*i; 

K(ll) = 0.266666667; 

Lh(ll) = -0.3798-5.1084*i; 
tolpha(ll) = 0.375-3.75*i; 

Inlpha(ll) = -20.0337-3.6847*i; 

K(12) = 0.24; 

Iii(12) = -0.552-5.8242*i; 
tfilp)n(12) = 0.375-4.16667*1; 

Lalplia(12) = -25.319-3.5256*i; 

K(13) = 0.2; 

Ui (13) = -0.886-7.2760*i; 

Maltlia(13) = 0.375-5.0*i; 

Lalplva(13) = -37.766-2.8460*1; 

Mi = 0.5; 

% Ncn-dimaisional coefficients 

q = 0; 
vlwn = 1.6; 
tin = 4; 
a = -0.3; 
xaljlia = 0.2; 
ralfiia2 = 0.25; 

% Flutter Determinant 

% Search for real root using bisection method 
% Outside for loop conducts the freqiuency svreep 


for m = 1:13 

X(l) = 0.2; 

X(3) = 2; 

X(2) = (X(3)-X(l))/2; 
for its = 1:50 

A = nit*(l-X*(vAiwa^2)*(li-i*g))+Iii(m); 

B = uu*xalplv=i>Lalpha(m)-tli(m)*(0.5+a)} 

D = in i *xalpha+Mi-Lh (m) * (0.5+a); 

E=tni*ralE)ha2* (l-X*(l+i*g) )-0.5*(0.5+a)4felpha(m)-Laipha(m) *(0.5+a)+Ih(m) *(0.5+a 

* 2 ; 

DET = A.*E - B*D; 

R a real (DET); 
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if R(1)*R(2) < 0 
X{3) = X(2)J 
else 

X(l) = X(2)} 
aid 

X(2) = (X(3)4X(l))/2; 

ati 

Rroot (m) = X(2); 

Rerror(m) = X(3) - X(l)‘ 

% Search for imaginary root 

X(l) = 0.2; 

X(3) = 2; 

X(2) = (X(3)~X(1) )/2{ 
for its = 1:50 


A = nu* (1-X* (whwa /v 2) * (l+i*g)) +Lh(m); 

B = nu*xalpha+Lalpha (m) -Lh(m) 4 (0.5+a) * 

D a nu*xalpha+ffi-Lh(mj * (0.54a) i 

E=*u*ralpha2* (1-X* (Ui*g)) -0.5* (0.54a) 4klpha(m) -Uifcha(in) * (0* 54a)+l!i(in) * (0.54a 
~ 2 ; 

GET = A.*E - B*D; 

I = inBg(DET) ; 
if 1(1)*1(2) < 0 
X(3) = X(2) ,* 
else 


X(l) = X(2); 
ax! 

X(2) = (X(3)4X(1))/2) 
ari 

Iroot (m) = X(2 )} 

Ienor(m) = X(3) - X(l),* 

end 

Rrrot 

Iroot 

cliff = Rroot - Iroot; 

% Plot of real and imaginary roots 


4 


%plot(l ./ K, Rroot,1 ./ K, Iroot) 
%grid 

fedabelCl/k*) 
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%ylabel('(X)*) 

"tpause 

% Search of roots for critical frequency 


n = length (Rroot),* 
for m=l:n 

if diff(m) > 0 
krf=K(m); 
break 
end 
ad 

% Calculation of the flutter coefficent that lies between 
* m and m-1 reduced frequency. 

IovK(l) = 1/K(m-1); 

Xreal(l) = sqrt (Rroot (m-1)) ,* 

Xixnag(l) = gqrt(Iroot (m-1)) } 

IovK(2) = 1/K(m); 

Xreal(2) = gqrt (Rroot (m))‘ 

Ximag(2) = gqrt (Iroot (m)) J 

xl = IcwK(l) :0.0001‘IovR(2) ,• 
n = length(xl); 

yl = l/(IovK(2) - lovK(l)) * ((IovK(2) - xl)*Xredi(l) + (xl - IovK(l))*Xreal(2))‘ 
y2 = l/(IovK(2) - IovK(l))*({IovK(2) - xi)*Ximag(l) + (xl - Iovk(l))*XijiBg(2)); 

plot (xl,yl,xl,y2) 
grid 

xlabel('1/k') 
ylabel('gqrt(X)') 


for m = l:n 

if y2(m)-yl(tn) < 0, 
Iok = xl(m) 
Xroot = y2(m) 
break 
ad 
ad 

kcrit * 1/Iok 
coef = Iok / Xroot 


* 
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about LE; Yshift=50 
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Nondimensional Time 




















Nondimensional Time 























Unstable Airfoil - Kp=0.115; Pivot about LE; Yshift=50 















Kp=0.115; Pivot about LE; Yshift=0.75 
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Nondimensional Time 

































Wake Trace Kp = 1.0 Scale = 0.5 





















Vortex 












Wake Trace Kp(l) = 1.0 Kp(2) = 2.0 Scale = 0.5 



Vortex Position 
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